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Abstract. We present an efficient algorithm for simulating the time 
evolution due to a sparse Hamiltonian. In terms of the maximum degree 
d and dimension N of the space on which the Hamiltonian H acts for time 
t, this algorithm uses (d 2 (d + log* N) \\Ht\\) 1+a< - 1 ' queries. This improves 
the complexity of the sparse Hamiltonian simulation algorithm of Berry, 
Ahokas, Cleve, and Sanders, which scales like (d 4 (log*iV) \\Ht\\) 1+0<1) . 
To achieve this, we decompose a general sparse Hamiltonian into a small 
sum of Hamiltonians whose graphs of non-zero entries have the property 
that every connected component is a star, and efficiently simulate each 
of these pieces. 



1 Introduction 

Quantum simulation of Hamiltonian dynamics is a well-studied prob- 
lem [1~3] and is one of the main motivations for building a quantum 
computer. Since the best known classical algorithms for simulating quan- 
tum dynamics are inefficient, Feynman suggested that computers that are 
inherently quantum might be better at simulating quantum systems [4]. 
Besides simulating physics, Hamiltonian simulation has algorithmic ap- 
plications, such as adiabatic optimization [5], unstructured search [6], and 
the implementation of continuous-time quantum walks [7,8]. 

The input to the Hamiltonian simulation problem is a Hamiltonian 
H and a time t; the problem is to implement the unitary operator e~ lHt 
approximately. We say that a Hamiltonian acting on an iV-dimensional 
quantum system can be simulated efficiently if there is a quantum circuit 
using poly (log N, t,l/e) one- and two-qubit gates that approximates (with 
error at most e) the evolution according to H for time t. Since the time 
evolution depends on the product Ht, the size of the circuit should also 
be bounded by a polynomial in some quantity measuring the size of H. 
When H is sparse, most of its matrix norms have comparable values, so 
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the complexity of simulating H is not very sensitive to how its size is 
quantified. It is conventional to require that the scaling be polynomial in 
|| #||, the spectral norm of H. 

Lloyd presented a method for simulating quantum systems that can be 
described by a sum of local Hamiltonians [1]. A Hamiltonian is called local 
if it acts non-trivially on at most a fixed number of qubits, independent 
of the size of the system. 

This was later generalized by Aharonov and Ta-Shma [2] to the case 
of sparse (and efficiently row-computable) Hamiltonians. A Hamiltonian 
is sparse if it has at most poly(log N) nonzero entries in any row. It is 
efficiently row-computable if there is an efficient procedure to determine 
the location and matrix elements of the nonzero entries in each row. 

The complexity of this simulation was improved by Childs [9] and 
further improved by Berry, Ahokas, Cleve and Sanders [3]. Their algo- 
rithm has query complexity (<i 4 (log* N) \\Ht\\) 1+0 W, where d is the max- 
imum degree of the graph of the Hamiltonian H. These algorithms de- 
compose the Hamiltonian into a sum of Hamiltonians, each of which is 
easy to simulate. In this paper, we present a different method of de- 
composing the Hamiltonian, giving an algorithm with query complexity 

(d 2 (d + log*7V)ll#*ll) 1+o(1) - 

Note that the simulation of Ref. [3] has also been improved using a 
completely different approach [10, 11]. That algorithm is more efficient in 
terms of all parameters except the error e, on which its dependence is 
considerably worse. The algorithm we present here maintains the same 
dependence on e as in Ref. [3], providing the best known method for 
high-precision simulation of sparse Hamiltonians. 

2 Hamiltonians and graphs 

A Hamiltonian H acting on n qubits is a 2™ x 2 n Hermitian matrix. It 
can also be thought of as the weighted adjacency matrix of a graph on 
2 n vertices, where the weights are complex numbers and the weight of an 
edge from u to v is the complex conjugate of the weight of the edge from 
v to u. We call the undirected graph formed by connecting two vertices 
if and only if the edge between them has nonzero weight the graph of the 
Hamiltonian. 

A Hamiltonian is said to be d-sparse if it has at most d nonzero entries 
in each row (i.e., the maximum degree of its graph is d). We often associate 
properties of the graph of a Hamiltonian with the Hamiltonian itself. For 



2 



instance, we might say "H is a forest," meaning that the graph of H is a 
forest. 

A star graph is a tree in which one vertex (called the center) is con- 
nected to all the other vertices and there are no other edges. In other 
words, it is a complete bipartite graph K\ )T . We call a forest in which 
each tree is a star graph a galaxy. 

A directed graph is a directed forest (directed tree) if its undirected 
graph is a forest (tree) . A directed tree is an arborescence if it has a unique 
root v such that all edges point away from v. Alternately, there is exactly 
one directed path from v to any other vertex u. In an arborescence, the 
edges are always directed from the parent to the child. A directed forest 
in which each tree is an arborescence is called a forest of arborescences. 

We use several matrix norms in our analysis. These include the spec- 
tral norm, \\H\\ := max||j,|| =1 \\Hv ||; the maximum entry norm, max(H) : = 
maxjj \Hij\\ and the maximum column norm, mcn(if) := maxj ||/fej||, 
where ej is the j th column of the identity matrix. 

3 Problem description and previous results 

The problem is to approximately implement the unitary e~ lHt for a d- 
sparse and efficiently row-computable iV-dimensional Hamiltonian H for 
time t. As input, we are given black-box access to H, and the values of d, 
t, and N. Since the Hamiltonian is sparse and efficiently row-computable, 
there is a convenient black-box formulation of the problem that abstracts 
away the details of computing matrix entries and locations. The Hamil- 
tonian is provided as a black-box function /, which accepts a row index 
and an integer i £ {1,2, ... ,d} and outputs the column index and ma- 
trix element corresponding to the i th nonzero entry in that row, if one 
exists. More precisely, if the nonzero elements in row x are y\, y2, ■ ■ ■ , Ud x , 
where d x < d is the degree of x, then f(x,i) = (yi,H Xtyi ) for i < d x and 
f(x, i) = (x, 0) for i > d x . This black box can be implemented efficiently if 
the Hamiltonian to be simulated is sparse and efficiently row-computable. 

For each row x, we allow the order in which the yi are given by the 
oracle to be arbitrary (but fixed). We do not assume that there is a 
convenient ordering, such as the increasing order of labels. To use the 
black box in a quantum circuit, we define an equivalent unitary matrix 
Uf which performs the operation Uf \x, i,0) = \x, i, f(x, i)). 

Let us denote the minimum number of queries to Uf required to ap- 
proximately simulate e~ lHt (up to error e, as quantified by the trace dis- 
tance) by Q(H,t). A common approach to this problem breaks it into 
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two subproblems, which we call the Hamiltonian decomposition prob- 
lem and the Hamiltonian recombination problem. First the Hamiltonian 
is decomposed into a sum of easy-to-simulate Hamiltonians; then these 
Hamiltonians are simulated for short times in a specific manner so that 
the overall simulation is approximately the same as that of H. 

Since we will also follow the decomposition-recombination strategy, 
we review this approach as applied in Ref . [3] . The given Hamiltonian H 
is decomposed into a sum of m Hamiltonians, H = Y%Li Hj. Let Q{Hj) 
denote the number of queries required to simulate Hj for time t' given 
black-box access to H. In general, the number of queries required might 
depend on t', but in the simulations used here Q{Hj) is independent of t'. 
Note that Q(Hj) includes the number of queries required to decompose 
H into Hj as well as to simulate Hj. In Ref. [3], the Hamiltonians Hj 
are 1-sparse, and their decomposition uses 0(log* N) queries 4 to a black 
box for H. Since a 1-sparse Hamiltonian can be simulated with 2 queries 
given an oracle for the 1-sparse Hamiltonian [7,9], Q(Hj) = 0(log* N). 
More precisely, 

Theorem 1 (Hamiltonian edge decomposition [3]). If H is an N x 

N Hamiltonian with maximum degree d, then there exists a decomposition 
H = Y^jL^Hj, where each Hj is 1-sparse, such that m = 6d 2 and each 
query to any Hj can be simulated by making Q(Hj) = 0(log* N) queries 
to H. 

These Hamiltonians are then recombined using the Lie-Trotter for- 
mula, which expresses the time evolution due to H as a product of time 
evolutions due to the Hj . The unitary e~~ lHt is approximated by a product 
of exponentials e~ lHjt , such that the maximum error in the final state 
does not exceed e. We want to upper bound the number of exponentials 
required, iV CX p. Reference [3] proves the following. 

Theorem 2 (Hamiltonian recombination [3]). Let k be any positive 
integer. If H = Y^jLi Hj is a Hamiltonian to be simulated for time t by 
a product of exponentials e~ lH ^' , and the permissible error (in terms of 
trace distance) is bounded by e < 1 < 2m5 fc_1 \\H\\ t, then the number of 
exponentials required, -/V exp , is bounded by 
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Using the upper bound on the number of exponentials and the num- 
ber of queries needed to simulate any exponential, the total number 
of queries needed to simulate the Hamiltonian H satisfies Q(H,t) < 
iV CX p x max,- Q(Hj). With Q(Hj) = <3(log* N) and m = 6d 2 , we get 

Q(H, t) = ^5 2fc d 4 (log* N) \\H\\ t (^Mll^j j . ( 2 ) 

We see that Q(H, t) is almost linear in t, which is almost optimal due 
to a no-fast-forwarding theorem [3]. However, the dependence on d is not 
optimal. In the present paper we improve the dependence on d without 
affecting the other terms. The dependence on d has been improved in 
other approaches, but only at the expense of a worse dependence on the 
error e [10, 11]. 

We propose a new algorithm for solving the Hamiltonian decomposi- 
tion problem. This strategy breaks up the Hamiltonian into only m = 6d 
parts, but increases Q(Hj) to 0{d + log* N), improving the overall de- 
pendence on d and N. 



4 Hamiltonian decomposition 

The Hamiltonian decomposition problem is the problem of decomposing 
a Hamiltonian H into a sum of m Hamiltonians Hj such that given a label 
1 < j < m and a time t', the unitary e~ lH ^ can be efficiently simulated. 

We solve this problem by decomposing the Hamiltonian into m = 6d 
galaxies. To achieve this, we first decompose the given graph into d forests 
using the forest decomposition technique of Paneconesi and Rizzi [12]. The 
idea is to assign one of at most d colors to each edge of the graph (not 
necessarily a proper edge coloring) such that the edges of any particular 
color form a forest. Not only is this decomposition possible, but it has 
some special properties that are required later in Lemma 2. 

Lemma 1 (Forest decomposition). For any Hamiltonian H of maxi- 
mum degree d, there exists a decomposition H = J2t=i H c and an assign- 
ment of directions to the edges such that each H c is a forest of arbores- 
cences. Furthermore, given a color c and a vertex v, we can determine 
v 's parent in H c with one query (or determine that it is a root) and with 
O(d) queries we can determine the list of edges in H c incident on v. 

Proof We first describe a procedure that assigns a color c to each edge. 
H c then consists of all edges colored c. To color the edges, every vertex 
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proposes a color for each edge incident on it using the oracle in the fol- 
lowing way: if f(x, i) = (y, H X;y ), then x proposes color i for the edge xy. 
Similarly, y proposes a color for the edge xy. The edge is now colored 
using the proposal of the vertex with higher label (i.e., if x > y then the 
edge xy is colored with x's proposal). This coloring uses d colors, which 
is optimal up to constants since a d-sparse graph can have dn/2 edges, 
but forests have at most n — 1 edges. 

Now we assign directions to the edges and show that each H c has 
no cycles, which shows that each H c is a directed forest. The edge xy 
is directed from x to y if x < y. This choice of directions results in a 
directed acyclic graph, which has no directed cycles. To rule out non- 
directed cycles, we note that any such cycle must contain a vertex v for 
which both the edges of the cycle point toward v. This means the label of 
v is greater than that of its two neighbors. Thus the color of these edges 
was decided by v, which cannot happen since vertices propose different 
colors for different edges. 

To show that each tree in H c is an arborescence, we show that it has 
a unique root. Observe that a directed tree with more than one root must 
have a vertex with more than one parent. This again leads to the situation 
where a vertex has two incoming edges of the same color, which is not 
possible since these edges are colored by this vertex's proposal. 

To show that the parent of a vertex can be determined with one query, 
note that if p v is the parent of vertex v in H c , then the edge from p v to 
v must be directed toward v. Thus the color of this edge is decided by v. 
If this edge is in H c , it is colored c. So if v has a parent, it must be the 
c th neighbor of v. With one query to the oracle, we can determine the c th 
neighbor of v. If there is no such neighbor, this vertex has no parent and 
is a root in H c . Otherwise the output contains the label of the parent. 

Finally, we show how to determine the list of edges in H c incident on 
x with 0{d) queries. First we query the oracle at most d times to get the 
labels of all the neighbors of x. For a neighbor y where y < x, the edge 
between x and y is colored by c only if y is x's parent in H c . Thus we 
can discard all edges xy where y < x but y is not the parent of x. When 
y > x, an edge between x and y is colored c only if x is y's parent in 
H c , and it takes one query to verify this for each y. Thus with at most d 
additional queries we can determine if all such edges are colored c. □ 

This lemma shows how to decompose a Hamiltonian into directed 
forests. Let T be the Hamiltonian of such a forest. We will decompose 
T into a sum of 6 galaxies, T = T\ + T2 + • • • + Tq. This is achieved by 
using an extension of the "deterministic coin tossing" protocol of Cole and 
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Vishkin [13] by Goldberg, Plotkin and Shannon [14]. Their protocol gives 
a proper vertex coloring of an arborescence using only 6 colors making 
0(log* N) queries. Vertex coloring a directed forest of arborescences gives 
a galaxy decomposition of the forest, since all the edges that point to 
vertices of a particular color form a galaxy. 

Lemma 2 (Vertex coloring a forest). IfT is a forest of arborescences, 
and the parent of a vertex can be determined with one query to an oracle 
for T, then there exists a proper vertex coloring of T using 6 colors, 
such that the color of any vertex can be determined by making 0(log* N) 
queries. 

Proof. We first describe the vertex-coloring procedure for the forest. A 
simple observation is that we already possess a vertex coloring of the 
forest: the labels of the vertices. This is a trivial proper vertex coloring 
using N colors. Now we use a procedure that decreases the number of 
colors used by a logarithmic factor. Then we can run several rounds of 
this procedure to decrease the number of colors down to 6. Let Cj(x) be 
the color assigned to vertex x at the beginning of the j th round of the 
procedure. At the beginning of the first round, we have c±(x) = x. 

Let x be a vertex with parent p x . Assume that we started with a 
proper vertex coloring at the beginning of round j. Since we have a proper 
coloring, Cj(x) / Cj{p x ). Let k be the index of the first bit at which x 
and p x differ, and let b be the value of the k th bit of x. The new color for 
vertex x is the concatenation of k and b, denoted (k, b). If x is the root, 
we take k = 0. We claim that if each vertex performs this procedure, the 
result is a proper vertex coloring. 

For a contradiction, suppose there are two adjacent vertices that have 
been assigned the same color in round j. Without loss of generality, one 
of them is the parent of the other, so let them be y and its parent p y . 
Since we started with a proper coloring at the beginning of round j, 
Cj(y) / Cjipy), but now c j+1 (y) = c j+1 (p y ). Let c j+1 (y) = (k,b) where 
by definition k is the bit at which Cj(y) and Cj{p y ) differ, and b is the 
value of the k th bit of Cj(y). Since Cj + \{p y ) also equals (k,b), the k th bit 
of Cj{p y ) is b. But Cj(y) and Cj{p y ) are supposed to differ at the k th bit, 
so this is a contradiction. Therefore the coloring procedure is valid. 

It remains to show that if the colors of the vertices are updated in this 
way, we reduce the number of colors to 6 in 0(log* N) rounds. If Lj is the 
number of bits used to represent colors at the beginning of round j, then 
Lj+i = [log(Lj)] +1. Initially, L\ = [log AT] . This recurrence relation can 
be solved to yield Lj < 3 when j = log* N [14]. Further rounds cannot 
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decrease Lj below 3, since Lj + \ = Lj when Lj = 3. A length of 3 bits 
allows the use of only 8 colors. Now we run the procedure once more. 
Since there are 3 possible values for k, and 2 for b, there are at most 6 
different colors. The total number of rounds is now log* N + 1. 

To show that the color of a vertex can be determined with 0(log* N) 
queries, we note that the color of vertex x at the end of the first round 
depends solely on x and p x . In general, the color of vertex x at the end 
of j rounds depends only on its first j ancestors. To determine x's color 
after log* N + 1 rounds, we need the labels of its log* N + 1 ancestors, 
which can be found with log* N + 1 queries, since the parent of a vertex 
can be found with one query. □ 

We have shown that a Hamiltonian can be decomposed into d forests 
of arborescences, each of which can be vertex-colored with 6 colors. If we 
consider all the edges of one of the d forests that point to a vertex of a 
particular color, this graph is a galaxy. So this decomposes the original 
Hamiltonian into 6d galaxies. For this particular decomposition of the 
Hamiltonian to be useful, we need to show that galaxies can be simulated 
easily. 

Theorem 3 (Galaxy simulation). If Hj is a Hamiltonian whose graph 
is a galaxy of maximum degree d, and the oracle can identify which vertices 
are centers of stars, then the unitary operator e~ lH ^ can be simulated 
using 0(d) calls to an oracle for Hj. 

Proof. The key idea is that given a vertex v, we can learn everything 
about the star to which v belongs in 0(d) queries. If v is the center of the 
star, the oracle identifies it as the center, so we can query all its neighbors 
to learn everything about the star with at most d queries. If v is not the 
center, we can determine the center, which is the only neighbor of v, with 
only one query, and then learn the rest of the star with at most d queries. 

Let R(x) denote all the information about the star to which x belongs: 
the label of the center, the labels of the other vertices in some fixed order, 
and the weights of all the edges. It is essential that R(x) depend only the 
star and not the particular vertex x chosen from the star, so that if x 
and y belong to the same star then R(x) = R(y). Since we know that 
R(x) can be computed with 0(d) queries, we can implement the unitary 
U given by U\x, 0) = \x,R(x)} with 0(d) queries. 

The Hamiltonian we are trying to simulate, Hj, is a galaxy. Thus, if 
c is the center of a star, and its neighbors are yi with edge weights wi, 
then Hj |c) = J2i w i \Vi)- H x is not the center of a star, and the edge 
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between x and the center c has weight w x , then Hj \x) = w x |c). Let K be 
a Hamiltonian which is similar to Hj, but acts on the input state \x, R(x)} 
instead of \x). That is, K \c, R(c)} = J2 i w i \Vii R{Ui)) when c is the center, 
and K \x,R(x)) = w x \c,R(c)) otherwise. Note that although the second 
register looks different, it is unaffected by K since R(x) depends only 
on the star and not the vertex. Combining K with the unitary U above, 
we see that Hj = U^KU. In words, U first computes all the information 
about the star in another register, K performs the required Hamiltonian, 
and the U' uncomputes the second register, which was unaffected by K. 

This simulation is efficient since K can be simulated efficiently. More 
importantly for our purposes, K requires no queries to implement, since 
all the information about the star is already present in the second register. 
Thus the operation Hj = U^KU requires only as many queries as U and 

require, which is 0(d). □ 

Combining Lemma 1, Lemma 2, and Theorem 3 gives our Hamiltonian 
decomposition theorem. 

Theorem 4 (Hamiltonian star decomposition). There exists a de- 
composition H = X)j=i Rj> where each Hj is a galaxy, such that m = 6d 
and each galaxy Hj can be simulated for time t' using Q(Hj) = 0(d + 
log* N) queries to an oracle for H. 

Proof. From Lemma 1, Lemma 2, and Theorem 3, we know that the 
claimed decomposition is possible. It remains to show that any Hj can 
be simulated for time t' using 0(d + log* N) queries. 

To show this, let us implement Hj on the basis state \x). If the imple- 
mentation is correct on all basis states, it is correct for all input states by 
linearity. We are given 1 < c < d and 1 < t < 6, which together form the 
index j. We want to simulate the galaxy formed by edges in H c directed 
toward vertices colored t by the vertex coloring algorithm of Lemma 2. 

From the proof of Theorem 3, it is clear that if we can compute R(x), 
then we can implement U, and thereby simulate the desired Hamiltonian. 
R(x) contains all the information about the star to which x belongs. Using 
the result of Lemma 1, we can determine the list of x's neighbors in H c 
using 0(d) queries. By the result of Lemma 2, with 0(log* N) queries we 
can determine x's color according to the vertex coloring algorithm. 

If x's color is not t, then x must be the center of a (possibly empty) 
star in Hj. The only edges in this star point toward vertices of color 
t, so we compute the colors of all the children of x in H c . These can 
be computed using only the labels of their log* N + 1 nearest ancestors, 
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which are all common ancestors. Thus we can compute the colors of all 
of x's children using 0(log* N) queries in total. Now we know the star 
around x, and thus R(x), using 0(d + log* N) queries. 

If x's color is t, then x's parent is the center of star. The parent of 
x , p x , can be determined with one query. Since x and p x are in the same 
star, R(x) = R(p x )- Since p x is the center of a star, we can compute R(p x ) 
as described above; thus we can also compute R(x). 

We have shown that for any x, we can compute R(x) with 0(d + 
log* N) queries. Thus the unitary U in the proof of Theorem 3 can be 
simulated with 0{d + log* N) queries. By Theorem 3, this means we can 
implement Hj with 0(d + log* N) queries, as claimed. □ 

Now we can use our Hamiltonian decomposition theorem with the 
Hamiltonian recombination theorem (Theorem 2). Since we have Q(Hj) = 
0{d + log* N) from Theorem 4 and m = 6d from Lemma 1 and Lemma 2, 
we get our final result using Q(H,t) < N cxp x maxj Q(Hj): 

Q(H, t) = (^ 2k d\d + log* N) \\H\\ t (^EU^j 1/2k ^j . ( 3 ) 

When compared with the query complexity of (2), we see that this 
improves the scaling with d. Furthermore, when d = I?(log*iV), which 
is likely to be the case when d is not constant, (3) has no log* N term: 
the scaling (in terms of d and N) is (<i 3 ) 1+0 ( 1 ), as compared to (2) which 
scales like (d 4 log* AT) 1+ °( 1 ). 

5 Remarks and conclusion 

So far, we have measured the size of H using the spectral norm \\H\\. 
However, if we express the simulation complexity in terms of a different 
norm, then both Theorem 2 and equation (3) can be improved to give 
slightly better bounds. 

In the proof of Theorem 2, \\H\\ is used as a simple upper bound for 
maxj However, omitting this step gives a slightly stronger version 

of Theorem 2 with \\H\\ replaced by maxj \\Hj\\. For a 1-sparse Hamil- 
tonian, \\Hj\\ = max(Hj) < max(if) [15], so ||-£f|| can be replaced by 
max(H) in (2). However, this also leads to an improvement of (3). When 
Hj is a galaxy, \\Hj\\ = mcn(Hj) [15], and since Hj is entry-wise upper 
bounded by H, mcn(Hj) < mcn(H). Thus \\H\\ can be replaced with 
mcn(iJ) in (3). To directly compare the two simulations, we can ap- 
ply the bound men (H) < \fd max(iJ) [15] to express both query com- 
plexities in terms of msx(H). In these terms, we still find that star 
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decomposition improves over edge coloring: our algorithm uses at most 
{d 2 - 5 (d+log* N) max(#i)) 1+ °W queries, whereas the algorithm of Ref. [3] 
scales like (d 4 (log* N) max(#t)) 1+ °W . 

In conclusion, we have described a Hamiltonian decomposition tech- 
nique that reduces the query complexity of simulating sparse Hamiltoni- 
ans. By the degree-dependent lower bounds established in Ref. [15], we 
know that query complexities scaling like \\H\\ or \fdmax.(H) cannot be 
achieved. It would be interesting to see if the Hamiltonian decomposition- 
recombination framework can be used to further reduce the dependence 
on d while keeping a similar dependence on the error e, or to establish 
stronger limitations on the simulation of sparse Hamiltonians taking error 
dependence into account. 
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